---
title: "Implicitization for the spherical trochoid"
author: "Stéphane Laurent"
date: '2023-09-26'
tags: R, rgl, maths, geometry
rbloggers: yes
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    highlight: kate
    keep_md: no
highlighter: pandoc-solarized
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE)
```

If you didn't read or if you don't remember my  post 
[Using implicitization to split a ball](https://laustep.github.io/stlahblog/posts/giacR02.html),
take a look at it before reading this one. 
To sum up, I took a spherical curve, namely a *satellite curve*, and I 
derived a surface whose intersection with the sphere is the satellite curve. 
I did that thanks to the so-called implicitization process, based on 
Gröbner bases, and which can be performed with the **giacR** package.

What I like with this process is that it's almost impossible to guess the 
shape of the derived surface when a spherical curve is given. For example, 
I still played with the satellite curves and using this process for 
different values of the parameters, I obtained this funny animated surface:

![](./figures/doubleTrumpet.gif)

You can find the code which generates this animation in
[this gist](https://gist.github.com/stla/c48977956eea1cf1cd581c6a5eab7686).

Now we will see what happens when applying this process to a 
*spherical trochoid*, a very interesting curve.

The parametric equations of the spherical trochoid are
$$
\left\{\begin{aligned}
x(t) & = \bigl(qb-b\cos(\omega)+d\cos(\omega)\cos(qt)\bigr)\cos(t)+d\sin(t)\sin(qt) \\
y(t) & = \bigl(qb-b\cos(\omega)+d\cos(\omega)\cos(qt)\bigr)\sin(t)-d\cos(t)\sin(qt) \\
z(t) & = \sin(\omega)\bigl(b-d\cos(qt)\bigr)
\end{aligned}\right..
$$

As suggested by its name, the spherical trochoid is a spherical curve. 
Namely, it is a curve lying on the sphere with center $(0, 0, h)$ and radius 
$R$ given by 
$$
h = \frac{b-a\cos(\omega)}{\sin(\omega)}
\qquad \textrm{and} \qquad
R = \sqrt{a^2+h^2+d^2-b^2}
$$
where $a = qb$. 

Type *"spherical trochoid"* on Google or another web search engine, you will 
find some valuable information, including some animations showing the 
geometric construction of the spherical trochoid.

Let's write these equations in R with $q = 3$ as a function of $t$, and let's 
do a plot:

```{r defineTrochoid}
sphericalTrochoid <- function(t) {
  A <- cos(omega)
  B <- sin(omega)
  f <- 3*b - b*A + d*A*cos(3*t)
  cbind(
    f * cos(t) + d * sin(t)*sin(3*t),
    f * sin(t) - d * cos(t)*sin(3*t),
    B * (b - d*cos(3*t))
  )
}
# parameters
omega <- 1.4
b <- 5
d <- 10
# sampling
t_ <- seq(0, 2*pi, length.out = 360L)
stSample <- sphericalTrochoid(t_)
# we will plot the supporting sphere as well
a <- 3*b
h <- (b - a*cos(omega)) / sin(omega) # altitude of center 
R <- sqrt(a^2 + h^2 + d^2 - b^2)     # radius
```

```{r plotTrochoid, eval=FALSE}
library(rgl)
sphereMesh <- cgalMeshes::sphereMesh(z = h, r = R, iterations = 5L)
open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.7)
shade3d(sphereMesh, color = "yellow", alpha = 0.2, polygon_offset = 1)
lines3d(stSample, lwd = 3)
```
![](./figures/sphericalTrochoid.gif)

Note that $q=3$ corresponds to the number of loops.

Now, let's proceed as in the post 
[Using implicitization to split a ball](https://laustep.github.io/stlahblog/posts/giacR02.html): 
using the Gröbner implicitization implemented in the **giacR** package, we 
shall derive an equation of an implicit surface whose intersection with the 
sphere is the spherical trochoid. We provide the polynomial expressions of
$\cos(3t)$ and $\sin(3t)$ in function of $\cos(t)$ and $\sin(t)$, which can 
be obtained with **giacR** as well.

```{r giacR}
library(giacR)

equations <- paste(
  "x = (3*b - b*A + d*A*cos3t)*cost + d*sint*sin3t",
  "y = (3*b - b*A + d*A*cos3t)*sint - d*cost*sin3t",
  "z = B*(b - d*cos3t)", 
  sep = ", "
)
relations <- paste(
  "A^2 + B^2 = 1", 
  "cost^2 + sint^2 = 1",
  "cos3t = 4*cost^3 - 3*cost", 
  "sin3t = (4*cost^2 - 1)*sint", 
  sep = ", "
)
variables <- "cost, sint, cos3t, sin3t"
constants <- "A, B, b, d"

giac <- Giac$new()
results <- giac$implicitization(equations, relations, variables, constants)
giac$close()
length(results)
```

Thus, **giacR** provides $24$ polynomial expressions, each equal to $0$. I 
don't know what are all these expressions. The last one is nothing but 
`"A^2+B^2-1"`. I tried a couple of them, I found that the eleventh one 
provides what we're looking for, and then I retained this one. It is a long 
expression:

```{r giac11th}
( fxyz <- results[11L] )
```

This expression is equal to $0$, and this provides an isosurface whose 
intersection with the above sphere is the spherical trochoid. Let's plot it.

```{r isosurface, eval=FALSE}
library(rmarchingcubes)
A <- cos(omega)
B <- sin(omega)
f <- function(x, y, z) {
  eval(parse(text = fxyz))
}
n <- 300L
x_ <- y_ <- seq(-R*1.1, R*1.1, length.out = n)
z_ <- seq(h - R - 0.1, h + R + 0.1, length.out = n)
Grid <- expand.grid(X = x_, Y = y_, Z = z_)
voxel <- with(Grid, array(f(X, Y, Z), dim = c(n, n, n)))
surf <- contour3d(voxel, level = 0, x_, y_, z_)
mesh <- tmesh3d(
  vertices = t(surf$vertices),
  indices  = t(surf$triangles),
  normals  = surf$normals
)
# plot
open3d(windowRect = 50 + c(0, 0, 512, 512))
shade3d(mesh, color = "darkviolet")
```
![](./figures/sphericalTrochoid_mesh.gif)

We get our surface! We will plot it along with the sphere and the curve. But 
I prefer to round its corners before, by clipping it to a horizontal cylinder.

```{r isosurfaceClipping, eval=FALSE}
# clip the mesh to a cylinder
fn <- function(x, y, z) x*x + y*y
cmesh <- clipMesh3d(mesh, fn, bound = (R*1.1)^2, greater = FALSE)
# plot everything
open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.75)
shade3d(sphereMesh, color = "orange", alpha = 0.4, polygon_offset = 1)
shade3d(cmesh, color = "darkviolet", polygon_offset = 1)
lines3d(stSample, lwd = 3, color = "chartreuse")
```
![](./figures/sphericalTrochoid_plotall.gif)